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Abstract 



" " I ' A new numerical model of the nonlinear diffusive shock acceleration is presented. It is used for modeling of particle acceleration in 
^—1 supernova remnants. The model contains coupled spherically symmetric hydrodynamic equations and the transport equations for 
energetic protons, ions and electrons. The forward and reverse shocks are included in the consideration. The spectra of cosmic rays 
released into interstellar medium from a supernova remnant are determined. The role of the reverse shock in the production of CR 
ions and positrons is discussed. 
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1. Introduction 

The diffusive shock acceleration (DSA) process |[ll|2,[3l01 is 
considered as the principal mechanism for production of galac- 
,tic cosmic rays (CR) in supernova remnants (SNRs). A signifi- 
cant theoretical progress in the investigation of this mechanism 
Was achieved (see e.g. Malkov & Drury |5] for a review). How- 
ever only during the last decade the excellent results of X-ray 
and gamma-ray astronomy supplied the observational evidence 
of the presence of multi-TeV energetic particles in these objects 
,(see e.g. Aharonian et al. [[60). 

■ Two shocks are produced by super-sonically moving super- 
nova ejecta after an explosion. A forward shock propagates in 
the circumstellar medium while a reverse shock propagates in 
the gas of ejecta. Generally it is believed that some part of ther- 
mal particles is injected at the shock fronts into acceleration. 

In this paper we present a new numerical model of nonlinear 
diffusive shock acceleration. This model is a natural develop- 
ment of the existing models QiMi- The solution of spherically 
symmetric hydrodynamic equations is combined with the en- 
ergetic particle transport and acceleration at the forward and 
reverse shocks of a supernova remnant. Nonlinear response of 
energetic particles via their pressure gradient results in the self- 
regulation of acceleration efficiency. 

Our previous studies that used this model dealt with the CR 
spectra produced by SNRs [91 11^. The input of the reverse 
shock was not taken into account there while the acceleration 
by both shocks was considered for modeling of non-thermal 
electromagnetic emission from the SNR RX J 1 7 1 3 .7-3946 HI | . 

The paper is organized as follows. The short description of 
the model is given in Sect. 2. The results of modeling of evolu- 
tion of a supernova remnant in the interstellar medium are pre- 
sented in Sect. 3. Sect. 4 contains the discussion of our results. 
Our conclusions are given in the last Section. The numerical 
code is described in Appendix. 



2. Nonlinear model of diffusive shock acceleration 

Hydrodynamical equations for the gas density p(r, f), gas ve- 
locity u{r, t), gas pressure Pg{r, t), and the equation for isotropic 
part of the CR proton momentum distribution N(r, t, p) in the 
spherically symmetrical case are given by 
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Here Pc - 4-n J p^dpvpN 13 is the CR pressure, w(r, f) is the ad- 
vective velocity of CRs, jg is the adiabatic index of the gas, and 
D{r, t, p) is the CR diffusion coefficient. It was assumed that 
the diffusive streaming of CRs results in the generation of mag- 
netohydrodynamic (MHD) waves. CR particles are scattered 
by these waves. That is why the CR advective velocity w may 
differ from the gas velocity u. Damping of these waves results 
in an additional gas heating. It is described by the last term in 
Eq. (3). Two last terms in Eq. (4) correspond to the injection 
of thermal protons with momenta p - Pf, p - Pb and mass m 
at the fronts of the forward and reverse shocks at r = Rfit) and 
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r - Rb(t) respectivel}Q. The dimensionless parameters 77^ and 
7/* determine the injection efficiency. 

The equation for ions is similar to Eq. (4). For ions with the 
mass M = Am and the mass number A it is convenient to use 
the momentum per nucleon p and the normalization of the ion 
spectra A^, to the baryonic number density. Then the number 
density of ions «/ is n, - \nA-^ / p^dpNi. The ion pressure 
Pi - 4-n J p-^dpvpNi/3 is also taken into account in the cosmic 
ray pressure Pc- 

We shall neglect the pressure of energetic electrons. In other 
words they are treated as the test particles. The evolution of 
electron distribution is described by equation similar to Eq. 
(4) with additional terms describing synchrotron and inverse 
Compton (IC) losses. 

CR diffusion is determined by magnetic inhomogeneities. 
Strong streaming of accelerated particles changes medium 
properties in the shock vicinity. CR streaming instability re- 
sults in the high level of MHD turbulence ^ and even in the 
amplification of magnetic field in young SNRs fl5]. Due to 
this effect the maximum energy of accelerated particles may be 
higher in comparison with previous estimates of Lagage and 
Cesarsky |fl4|. 

According to the numerical modeling of this non-resonant 
instability, the magnetic field is amplified by the flux of run- 
away highest energy particles in the relatively broad region up- 
stream of the shock [13]. Magnetic energy density is a small 
fraction (~ 10"^) of the energy density of accelerated particles. 
This amplified almost isotropic magnetic field can be consid- 
ered as a large-scale magnetic field for lower energy particles 
which are concentrated in the narrow region upstream of the 
shock. Resonant streaming instability of these particles pro- 
duces MHD waves propagating in the direction opposite to the 
CR gradient. Strong nonlinear damping of these waves results 
in the gas heating (see the last term in Eq.(3)). CR gradient is 
negative upstream of the forward shock and MHD waves prop- 
agate in the positive direction. The situation changes down- 
stream of the forward shock where CR gradient is as a rule pos- 
itive and MHD waves propagate in the negative direction. This 
effect is mostly pronounced downstream of the forward shock 
of SNR because the magnetic field is additionally amplified by 
the shock compression and the Alfven velocity Va - B/ ^jAnp 
can be comparable with the gas velocity in the shock frame 
u' - Rf - u(Rf - 0,t). As for CR diffusion coefficient, it is 
probably close to the Bohm value Db = pvc/3qB, where q is 
the electric charge of particles. 

We apply a finite-difference method to solve Eqs (1-4) nu- 
merically upstream and downstream of the forward and reverse 
shock (see Appendix). A non-uniform numerical grid upstream 
of the shocks at r > Rf and r < Ri, allows to resolve small 
scales of hydrodynamical quantities appearing due to the pres- 
sure gradient of low-energy CRs. The gases compressed at the 
forward and reverse shocks are separated by a contact disconti- 
nuity aXr - Rc that is situated between the shocks. An explicit 



conservative TVD scheme 11511 for hydrodynamical equations 
(1-3) and uniform spatial grid were used between the shocks. 

Because of synchrotron losses the spacial scale of the high- 
energy electrons can be rather small downstream of the shocks. 
That is why we use a non-uniform numerical spacial grid for 
accelerated electrons downstream of the shocks. 

The magnetic field plays no dynamical role in the model. 
Since we have not perform the modeling of the amplification 
and transport of magnetic field here, its coordinate dependence 
should be specified for determination of cosmic ray diffusion 
and for calculation of synchrotron emission and losses. We 
shall assume below that the coordinate dependencies of the 
magnetic field and the gas density coincide upstream and down- 
stream of the forward shock: 
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Here po is the gas density and V^o = ^o/ V^^Po is the Alfven 
velocity of the circumstellar medium. The parameter Ma de- 
termines the value of the amplified magnetic field strength. For 
low shock velocities Rf < MaVao the magnetic field is not am- 
plified. 

The magnetic energy is about 3.5 percent of the dynamical 
pressure pqR^ according to estimates from the width of X-ray 
filaments in young SNRs ||l6ll . This number and characteristic 
compression ratio of a modified SNR shock <j - 6 correspond 
to Ma ~ 23. Since the plasma density p decreases towards 
the contact discontinuity downstream of the forward shock, the 
same is true for the magnetic field strength according to Eq. (5). 
This seems reasonable because of a possible magnetic dissipa- 
tion in this region. 

Situation is different downstream of the reverse shock at 
Rb < r < Re- The plasma flow is as a rule strongly influenced by 
the Rayleigh- Taylor instability that occurs in the vicinity of the 
contact discontinuity and results in the generation of MHD tur- 
bulence in this region. We shall assume that the magnetic field 
does not depend on radius downstream of the reverse shock 
while the dependence in the upstream region is described by 
the equation similar to Eq. (5): 



B{r, t) = V47rp, 
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We shall use indexes / and b for quantities corresponding to the forward 
and reverse (backward) shock respectively. 



Here r„ < Rb is the radius where the ejecta density has a min- 
imum and equals p„,. This radius r„, is generally close to the 
reverse shock radius Rb and is equal to it if the reverse shock is 
not modified by the cosmic ray pressure. 

CR advective velocity may differ from the gas velocity on 
the value of the radial component of the Alfven velocity VAr 
calculated in the isotropic random magnetic field: w = u + 
^aVa/ V3. Here the factor describes the possible deviation of 
the cosmic ray drift velocity from the gas velocity. The similar 
expression for the cosmic ray drift velocity is used upstream of 
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Figure 1: Dependencies on time of the forward shock radius Rf (thick solid 
line), the reverse shock radius R\, (thick dashed Hne), the forward shock velocity 
Vf (thin solid line) and the reverse shock velocity Y], (thin dashed line). The 
ratio of CR energy and energy of supernova explosion E^rlEsN (dotted line) is 
also shown. 
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Figure 2: Radial dependencies of the gas density (thick solid line), the gas 
velocity (dotted line), CR pressure (thick dashed line) and the gas pressure 
(dashed line) at f = 10^ yr. At this epoch the forward shock velocity is 3300 km 
s"', its radius is 6.5 pc, the reverse shock velocity is 1650 km s"', its radius is 
5. 1 pc, the magnetic field strength downstream of the forward shock is 160 [iG 
while the magnetic field strength downstream the reverse shock is 56 {iG. 



the reverse shock at r < Rb- We shall use values = 1 and = 
-1 upstream of the forward and reverse shocks respectively, 
where Alfven waves are generated by the cosmic ray streaming 
instability and propagate in the corresponding directions. The 
damping of these waves heats the gas upstream of the shocks 
(ill and limits the total compression ratios by a number close to 
6. This Alfven heating produce a more efficient limitation of the 
shock modification in comparison with the dynamical effects of 
the magnetic field considered by Caprioli et al. [ 1 8] which are 
neglected in our study. In the downstream region of the forward 
and reverse shock at Rb < r < R f we set =0 and therefore 
w = M in the major part of this paper except the consideration 
of the Alfven drift effects downstream of the shocks at the end 
of Sect.3 and in Fig. 10. 
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Here the parameter g > Q depends on the type of nonlinear 
wave damping which is relevant only for low velocity shocks 
Rf < MaVao when the magnetic field is not amplified. The 
parameter tjb describes the possible deviations of diffusion co- 
efficient from the Bohm value Db = vpc/3qB{r,t) for high- 
velocity shocks. Since the highest energy particles are scattered 
by small-scale magnetic fields, their diffusion is faster than the 
Bohm diffusion 1. 1 3.1 . The same is true for smaller energy par- 
ticles because they can be resonantly scattered only by a frac- 
tion of the magnetic spectrum. We shall use the value tjb - 2 
throughout the paper. 

We shall use the value of parameter g = 1 .5 below. It corre- 
sponds to the nonlinear wave damping of the weak turbulence 
theory. Note that a stronger Kolmogorov-type nonlinear damp- 
ing used by Ptuskin & Zirakashvili lll9ll for estimate of maxi- 
mum energy in SNRs corresponds to g -?>. 

In real situation the level of MHD turbulence drops with dis- 
tance upstream of the shock and diffusion may be even faster 
there. This is qualitatively taken into account by the exponen- 
tial factors in Eq. (7). The characteristic diffusive scale of high- 
est energy particles is a small fraction << 1 of the shock 
radius |13] and is determined by the generation and transport 
of MHD turbulence in the upstream region 1 2(1 ,2 1.] . The value 
^0 ~ In '(Z)c/£>j) is determined by ratio of diffusion coeffi- 
cient Dc in the circumstellar medium and diffusion coefficient 
Ds « Dc in the vicinity of the shock. The MHD turbulence 
is amplified exponentially in time before the shock arrival from 
the background level by cosmic ray streaming instability. The 
characteristic range of is 0.05 -^0.1 [il3il . We shall use ffie 
value ^0 - 0.05 below. 

It is believed that the supernova ejecta has some velocity dis- 
tribution P{V) just after the supernova explosion 12211 
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Here the index k characterizes the steep power-low part of this 
distribution. The radial distribution of ejecta density is de- 
scribed by ffie same expression with V - rjt. The character- 
istic ejecta velocity Yej can be expressed in terms of energy of 
supernova explosion E^f^ and ejecta mass M^j as 



10(^-5)£^jv 
3(;t-3)M, 



1/2 



(9) 



3. Numerical results 

Figures (l)-(9) illustrate the numerical results that are ob- 
tained for the SNR shock propagating in the medium with 
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Figure 3: Spectra of accelerated particles at f = 10^ yr. The proton spectrum 
at the forward shock (thick solid), ion spectrum at the reverse shock (thick 
dashed), electron spectrum at the forward shock (multiplied to the factor of 100, 
thin solid) and positron spectrum at the reverse shock (multiphed to the factor 
of 100, thin dashed) are shown. Spectrum of ions is shown as the function of 
momentum per nucleon and normalized to the baryonic number density. 



Figure 5: Spectra of accelerated particles at / = 10 yr. The spectrum of protons 
(thick solid line) and electrons (multiplied to the factor of 100, thin sohd hne) 
at the forward shock, ion spectrum (thick dashed line) and positron spectrum 
(multiplied to the factor of 100,thin dashed line) in the central part of the rem- 
nant are shown. Spectrum of ions is shown as the function of momentum per 
nucleon and normalized to baryonic number density. 




Figure 4: Radial dependencies of the gas density (thick solid line), the gas 
velocity (dotted line), CR pressure (thick dashed line) and the gas pressure 
(dashed line) at ( = 10** yr. At this epoch the forward shock velocity is 660 
km s"', its radius is 17 pc, magnetic field strength downstream of the forward 
shock is 40 fiG. 

a hydrogen number density hh - 0.1 cm"^, magnetic field 
strength Bo = 5 and temperature T - 10^ K. The frac- 
tion XHe - nHe/riH = 0.1 of hclium nuclei was assumed. The 
gas of ejecta does not contain hydrogen in the case considered. 
We use the ejecta mass M^j = 1 AMq, the energy of explosion 
EsN = 1-0 ■ 10^'' erg and the parameter of ejecta velocity dis- 
tribution k - 1 . The value of the parameter - 23 was 
assumed. 

The initial forward shock velocity is Vq = 2.9 ■ 10"* km 
s~'. The injection efficiency is taken to be independent on 
time rf - rj^ - 0.01, and the injection momenta are pf = 
2m(R f-u(R + 0, t)), ph - 2m(u(Rh - 0, f) - Rb)- Protons with a 
mass m are injected at the forward shock while ions with mass 
number A and charge number Z - A/2 are injected at the re- 
verse shock. This high injection efficiency results in the signif- 
icant shock modification already at early stages of SNR expan- 



sion while the thermal sub-shock compression ratio is close to 
2.5 during the simulation. This is in agreement with the radio- 
observations of young extragalactic SNRs ll23ll and with the 
modeling of collisionless shocks (2^. Similar values of the in- 
jection efficiency were found in hybrid modeling (see e.g. 1,25.] ) 
and at the Earth bow shock in the solar wind 02611 . 

As for the electron injection we assume a rather high injec- 
tion energy of electrons Ei„j = 100 MeV. This qualitatively 
corresponds to some models of suprathermal electron injec- 
tion. Partially ionized ions accelerated at the shocks up to 
relativistic energies can produce multi-MeV electrons in the 
upstream region in the course of photo-ionization by Galac- 
tic optical and infrared radiation MeV electrons and 
positrons are also present in the radioactive supernova ejecta 
while gamma-rays from ^^Co decay in ejecta produce ener- 
getic electrons via Compton scattering in the circumstellar 
medium |!28]. These energetic particles may be additionally 
pre-accelerated via stochastic acceleration in the turbulent up- 
stream regions of the shocks. 

Below we assume that electrons are injected at the forward 
shock with efficiency rf_ - 10 Rf /c while positrons are in- 
jected at the reverse shock with efliciency = 10"^. These 
numbers are expected for the injection mechanisms mentioned 
above (see Sect.4 for details). Since electrons are consid- 
ered as test particles our results can be easily rescaled for any 
other injection efficiency. The chosen injection rate at forward 
shock maintains the electron to proton ratio K^p of the order of 
K-p ~ 10 throughout the simulation while time-independent 
positron injection at the reverse shock results in positron to 
ion ratio K+i increase from K+i ~ 10 in the very beginning 
of SNR evolution up to K+j ~ 10"^ at several thousand years 
shortly before the disappearance of the reverse shock. 

In the real SNR the ions are also injected at the forward 
shock. We do not consider this process here in order to find 
the spectra of ions produced by the reverse shock. Because of 
the same reason we assumed the absence of hydrogen in the 
ejecta. Although this is the case for la/b/c and lib supernovae it 
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Figure 6: Radial dependencies of the gas density (thick solid line), the gas 
velocity (dotted line), CR pressure (thick dashed line) and the gas pressure 
(dashed line) at t = lO' yr At this epoch the forward shock velocity is 164 
km s" ' , its radius is 42 pc, magnetic field strength downstream of the forward 
shock is 23 fiG. 



is not true for IIP supernovae. It is expected that the spectra of 
ions injected at the forward and reverse shock are similar to the 
spectra of protons. The production of CR ions at the forward 
shock was recently considered by Caprioli et al. ll29ll . 

The dependencies on time of the shock radii Rf and Rh, the 
forward and reverse shock velocities Vy - Rf and Vb - Rh, 
CR energy Ecr/EsN are shown in Fig.l. The calculations were 
performed until f = 10^ yr, when the value of the forward shock 
velocity drops down to ^/ = 164 km s ' and the forward shock 
radius is Rf = 42 pc. 

At early times of SNR evolution the distance between reverse 
and forward shocks is only 10% of the remnant radius. This 
is less than 23% thickness for automodel Chevalier-Nadezhin 
solution with k - 7 [301 and should be attributed to a strong 
modification of both shocks by CR pressure. The reverse shock 
is strongly decelerated only when the forward shock sweeps 
the gas mass comparable with the ejecta mass at f > 100 yr and 
when the transition to the Sedov phase begins. 

Radial dependencies of physical quantities at f = 10^* yr are 
shown in Fig. 2. The contact discontinuity between the ejecta 
and the interstellar gas is at r - Rc = 5.6 pc. The reverse shock 
in the ejecta is situated at r - Rh = 5.1 pc. At the Sedov stage 
the reverse shock moves in the negative direction and reach the 
center after seven thousand years after the supernova explosion. 
We stop the calculations in the region r < Rh when the reverse 
shock radius Rh < 0.1/?/. 

It should be noted that our one-dimensional calculations 
cannot adequately describe the development of the Rayleigh- 
Taylor instability of the contact discontinuity. In the real situ- 
ation the supernova ejecta and the circumstellar gas are mixed 
by turbulent motions in this region (see e.g. MHD modehng of 
Jun & Norman |31]). 

Spectra of accelerated protons and electrons at f = 10-* yr 
are shown in Fig. 3. At this epoch maximum energy of protons 
accelerated in this SNR is about 100 TeV, while higher energy 
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Figure 7: Spectra of accelerated particles at / = 10^ yr The spectrum of protons 
(thick solid line) and electrons (multiplied to the factor of 100, thin sohd hne) 
at the forward shock, ion spectrum (thick dashed line) and positron spectrum 
(multiplied to the factor of 100, thin dashed line) in the central part of the rem- 
nant are shown. Spectrum of ions is shown as the function of momentum per 
nucleon and normalized to baryonic number density. 
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Figure 8: Spectra of particles produced in the supernova remnant during lO' yr. 
Spectrum of protons injected at the forward shock (thick solid line ), spectrum 
of electrons injected at the forward shock (thin solid line), spectrum of ions 
injected at the reverse shock (thick dashed line) and the spectrum of positrons 
injected at the reverse shock (thin dashed line) are shown. Spectrum of ions is 
shown as the function of momentum per nucleon and normalized to the bary- 
onic number density. 



particles have already left the remnant. 

Radial dependencies of physical quantities at later epoch t - 
IQ^ yr are shown in Fig.4. The reverse shock has reached the 
center of the remnant earlier A weak reflected shock is clearly 
visible at r = 5 pc. 

The spectra of particles ar f = 10"^ yr are shown in Fig. 5. 

Radial dependencies of physical quantities at the end of sim- 
ulation at t - 10^ yr are shown in Fig. 6. At this epoch the 
remnant is deeply in Sedov stage. The contact discontinuity is 
at r = 14 pc. 

The spectra of particles at f = 10^ yr are shown in Fig. 7. The 
shock modification is not strong because the Alfvenic Mach 
number is close to 6 and the corresponding Alfven heating up- 
stream of the forward shock results in the lower compression 
ratio and acceleration efficiency. That is why the spectra of par- 
ticles are steeper in comparison with ones at earlier epochs. 

The spectra of particles produced during the whole evolution 
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Figure 9: Dependencies on time of gamma-ray fluxes from SNR at distance 4 
kpc. We show the gamma-ray flux at 1 TeV produced via pion decay (thick 
solid line) and via IC process (thick dashed line), gamma-ray flux at 1 GeV 
produced via pion decay (thin solid line) and via IC process (thin dashed line). 
For the sake of simplicity only IC scattering of microwave background photons 
was considered. The evolution of synchrotron X-ray flux at 2 keV (thick dotted 
Une) and radio-flux at 1400 MHz f 1400 (dotted line) are also shown. 

of the remnant are shown in Fig. 8. They are obtained as the sum 
of the spectra integrated throughout simulation domain and of 
the time-integrated diffusive flux at the simulation boundary at 
r = 2Rf. At f = 10^ yr the maximum energy of currently ac- 
celerated particles drops down to 100 GeV because of nonlinear 
damping. Higher energy particles have already left the remnant. 
Note that stronger Kolmogorov-type damping with g = 3 will 
result even in lower energies of the order of 1 GeV. However 
we found that the spectra do not change in this case. 

We found that the maximum energy of CR protons is some- 
what less than 10'^ eV. It is almost an order of magnitude lower 
in comparison with the results of Ptuskin et al. 1 10] where more 
optimistic assumptions tjb - I and the spatially-uniform CR 
diffusion coefficient upstream of the forward shock were used 
(cf. Eq.(7)). 

Note that the synchrotron losses of run-away electrons and 
positrons were taken into account in our modeling. The cut- 
off energy of the leptonic spectra ~ 5 TeV shown in Fig. 8 is 
determined by the magnetic field strength Bo = 5yuG in the cir- 
cumstellar medium and by the remnant age f = 10^ yr 

Evolution of non-thermal emission produced in the SNR at 
distance 4 kpc is shown in Fig. 9. It is worth noting that a sig- 
nificant part of the IC emission is produced in the central region 
of the remnant at late epochs when the reverse shock have dis- 
appeared. This is because the magnetic field is rather weak in 
these regions. 

4. Discussion 

Although only about 5% of supernova energy is transferred 
to the particles accelerated at the reverse shock, they cannot be 
neglected. First of all the ejecta has absolutely different com- 
position in comparison with interstellar medium where the for- 
ward shock propagates. Since the solar abundance corresponds 
to 1% in the mass of heavy elements while the ejecta can con- 
tain up to 50% of heavy elements it is clear that the reverse 



Figure 10: Spectra of particles produced in the supernova remnant during 10^ 
yr in the model including the Alfven drift downstream of the shocks. Spectrum 
of protons injected at forward shock (thick solid line ), spectrum of electrons 
injected at the forward shock (thin solid line), spectrum of ions injected at the 
reverse shock (thick dashed line) and the spectrum of positrons injected at the 
reverse shock (thin dashed line) are shown. Spectrum of ions is shown as the 
function of momentum per nucleon and normalized to the baryonic number 
density. 

shock will dominate in the production of heavy high-energy CR 
nuclei. 

The relative input of the reverse shock at high energies is 
determined by the relative energetics of the reverse and forward 
shocks that is of the order of 1/10 for SNRs expanding in the 
uniform circumstellar medium. 

According to our results, approximately 70% of supernova 
energy is transferred to particles accelerated by forward shock. 
It is significantly higher than the estimate of 10-20 % needed 
to maintain CR density in the Galaxy if the supernova rate is 
1/30 yr One of possibilities to resolve this contradiction is the 
assumption that CRs are accelerated only at small part of the 
forward shock surface. This can be due to the dependence of 
the proton and ion injection on the shock obliqueness 1 32] . This 
effect is observed in SN 1006 and in the interplanetary medium. 

This effect does not influence strongly the ion injection at the 
reverse shock. It is expected that the random magnetic field is 
the main component of the field in the expanded ejecta. This 
is because the magnetic field of ejecta originates from the mag- 
netic field of the exploded star Thus the random magnetic field 
strength of the red super giant progenitor of IIP supernovae is of 
the order of 10^ G similar to the magnetic field strength in the 
Sun interior while the regular field is of the order of 1 G. After 
a homogenous expansion from the initial stellar radius 10'^ cm 
up to the radius 10'^ cm of a young SNR the frozen-in magnetic 
field drops down to 10"^ G. Although this value is significantly 
lower than the magnetic field in the interstellar medium it is 
strong enough for acceleration of particles up to 100 GeV in 
SNRs. A more realistic non-homogenous expansion will result 
in the stretching of the field in the radial direction. This can 
increase the magnetic field strength and ion injection efficiency 
at the reverse shock. Magnetic fields can be also amplified by 
non-resonant CR streaming instability suggested by Bell [12J. 
If so the relative contribution of the reverse shock to the over-all 
CR spectrum increases. 



6 



The effect also depends on the type of the supernova explo- 
sion. It is known that acceleration at the reverse shock occurs 



in Cas A SNR 11331 13411 . The progenitor of the core-collapse 



Cas A supernova had the radius of the order of 10'^ - 10'^^ cm. 
On the other hand white dwarfs that are progenitors of la super- 
nova explosions have small radii of the order of 10^ cm. Interior 
magnetic field 10^ G of the white dwarf will drop down to 10 
G after the homogenous expansion of the young SNR. Such a 
weak magnetic field can result in the ineffective DSA at the re- 
verse shock of la supernovae. This effect is probably observed 
in Tycho SNR (see Warren et al. llssll for details). 

Another possibility to suppress the CR production is related 
with the Alfven drift downstream of the forward shock |9]. It 
results in the steeper spectrum of CRs accelerated at the forward 
shock. On the other hand the Alfven drift downstream of the 
reverse shock may produce even an opposite effect because the 
CR gradient is positive in this region (see Fig. 2). As a result the 
input of the reverse shock will be significant at TeV energies. 

The over-all spectrum according to this model is shown in 
Fig. 10. The value of = -1 instead of = is used down- 
stream of the shocks. Because of the Alfven drift the positron 
spectrum at reverse shock is significantly harder than the elec- 
tron spectrum at the forward shock. 

It should be noted that the spectra of ions accelerated at the 
reverse shock are harder than the proton spectra at the forward 
shock (see Figs 3,8) in spite of the same level of the shock mod- 
ification for both shocks. This is because the shocks propagate 
in the media with different properties. When the reverse shock 
reach the flat part of the ejecta density distribution it propa- 
gates in the medium with decreasing in time density. That is 
why the number of freshly injected ions is low in comparison 
with higher energy particles accelerated earlier. It results in the 
spectral hardening. This effect is absent at the forward shock 
propagating in the medium with a constant density. The Alfven 
drift strengthens this effect (see Fig. 10). 

We adjust the electron (positron) injection to produce a suffi- 
cient number of electrons and positrons. The spectra shown in 
Figs 8 and 10 can explain CR Galactic electrons and positrons. 
The expected positron injection efficiency from the radioactive 
decay of ^^Ti is estimated as 77* ~ MnKAAM^j) ~ 10"* Izifor 
'^"^Ti mass of the order of ~ 10'^ Mq as observed in SNRs 13611 . 
We used this number in our simulation. As for the electron in- 
jection at the forward shock the relative number of energetic 
electrons from photo-ionization of accelerated single charged 
He ions is of the order of 77^ ~ xnej'^ in.^^(Pmax/fnc)R/ /c^ ~ 
lQ~^Rj,/c^. Here y ~ iHel^ph ~ 10 is the gamma-factor of 
single-charged He ion photo-ionized by galactic ultra-violet 
photons with energy €ph ~ 10 eV and Ine - 52 eV is the 
ionization potential of helium. This slightly overestimates the 
electron injection in young SNRs where the acceleration is fast 
enough, y is closer to y ~ 100 and the ionization is provided by 
eV optical photons. However we used this crude estimate ffiat 
is justified in old SNRs for electron injection in our simulations. 
This injection mechanism suggested by Morlino [27] produces 
one order of magnitude higher number of energetic electrons 
in comparison with the number 77^ ~ 10 ''7?^ of Compton 



scattered electrons energized by gamma-photons from ^*Co ra- 
dioactive decay in supernova ejecta ll28ll . Here R/^pc is ffie for- 
ward shock radius expressed in parsecs. 

The efficiency of electron (positron) acceleration at the re- 
verse shock in SNRs can be observationally checked in the 
gamma-ray band. Young and middle-age (t ~ 10'^ yr) SNRs can 
be bright in TeV energies (see Fig. 9). The electrons (positrons) 
accelerated before at the reverse shock do not strongly influ- 
enced by synchrotron losses in the weak magnetic field of the 
central part of the SNR and emit IC gamma-rays. 

5. Conclusion 

The main purpose of the present paper was the presentation 
of a new numerical code for the modeling of hydrodynamics 
and nonlinear shock acceleration in SNR. The code is described 
in details in five Appendixes. It develops the modeling of par- 
ticle acceleration by spherical shocks fulfilled earlier by other 
authors [7, 8]. Besides some important technical details, the 
main novel features of our code include the particle accelera- 
tion by two shocks - forward and reverse, the evolving with a 
SNR age magnetic field, which determines the cosmic ray trans- 
port, the account for the Alfven drift effects both upstream and 
downstream of the forward and reverse shocks. The first ver- 
sion of the model of acceleration and a number of important 
applications including the explanation of the overall spectrum 
of galactic cosmic rays ftlOil , and the modeling of particle ac- 
celeration and nonthermal radiation in SNR RX J 17 13.7-3946 
1 1 ill were developed upon our code refinements. 

The special attention in the astrophysical applications of the 
modeling discussed in the present work was focused on the par- 
ticle acceleration by the reverse SNR shock. It was shown that 
the reverse shock can give a non-negligible contribution to the 
production of CR ions and positrons as compared to the contri- 
bution of the forward shock. The spectra of particles acceler- 
ated at the reverse shock can be harder than the spectra at the 
forward shock, see Figures 8 and 10. It may offer a new inter- 
pretation of cosmic ray data [37] that suggests the presence of 
high energy primary positrons. The acceleration of positrons 
by the reverse shock moving through the ejecta material, which 
contains the positrons from "^^Ti radioactive decays, was pro- 
posed in [28]. It is also possible [38[ that the hard spectrum 
of accelerated nuclei and the low abundance of hydrogen in the 
material of supernova ejecta accelerated by the reverse shock 
may explain the difference in the energy spectra of hydrogen 
and helium in cosmic rays. We plan detailed study of all these 
effects in a future work. 

In the light of the problem of electron injection it is of inter- 
est that the models of suprathermal electron injection 11281 127 1 
reproduce the required amount of Galactic CR electrons and 
positrons if the leptons are pre-accelerated up to Ei„j ~ 100 
MeV in the upstream regions of supernova shocks. 

The work was supported by the Russian Foundation for Basic 
Research grant 10-02-001 10a. We thank the anonymous referee 
for a number of valuable comments. 



7 



Appendix A. Details of the numerical method 



We used the following change of variables upstream of the shocks: 



(A.1) 



and downstream of the shocks: 



R,-R„' Rh<r<Rc, 
ipt^Rc<r<Rf. 



(A.2) 



We shall use the dimensionless parameters i = \n{Rf/Ro) instead of time t and ^ = ln{p/mc) instead of momentum p. It 
is convenient to use the new variable n(^) - Ancp'^N(p) instead of momentum distribution N(p). For the relativistic momenta 
p » mc this variable corresponds to a partial energy density of cosmic rays. 

In these new variables the equations (l)-(4) have the following form upstream of the shocks: 



d 1 r, 5 T 1 Rf 



Rf d 



84 

d 



Vfdf 



d 1 1 d T I Rf o Rf 



(U - W)f—Pc + PgT-^^U 



(A.3) 



(A.4) 



(A.5) 



■< i^f dn d t dn dn ^ / /, 

Rluf — ^ = —RfbfD—-R]t,f{w - VfbO— + — 



dn 

— -An 



(A.6) 



Here indexes / and b are referred to the forward > 1) and reverse shock (0 < ^ < 1) respectively. 
The equations (l)-(4) have the following form downstream of the shocks: 



dt 



^ 2 f,b Rf 



(A.7) 



iRfl^P^fb d 



drj ' 



(A.8) 



d ,Rf 



,Rf d 



dt 



Vfdr] 



(A.9) 



,5n 1 



9 . Vf dn d r^D dn o th 

^' Rf dt dT]Afj,dT] dr] 3 



dn 

4n 

d^ 



drj 



2 

-r w. 



(A. 10) 



Here the speeds and are determined asu^ - u- Vdl - if)- tjVj and - u- Vdl +r]) + r/Vh respectively. The corresponding 
CR advective speeds and are determined as = w - VcC 1 - 77) - V/ and = w - VeC 1 + //) + r]Vb respectively. The quantities 
Af^b are the distances between the forward shock and the contact discontinuity A/ = Rf-R^ and between the reverse shock and the 
the contact discontinuity Ab = Rc - Rb- We introduced energy density of the gas e - ^ + Then the equations (A7)-(A9) are 
written in the conservative form that is convenient for the numerical method. The radius r in the last equations should be expressed 
using Eq. (A2): 



Rc + T](Rf-Rc), < 77 < 1, 
Rc + i](Rc - Rb), -1 <77<0. 



(A.11) 
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Appendix B. Solution of hydrodynamic equations in the upstream regions 



Since the flow upstream of the shocks is supersonic the hydrodynamic equations can be solved using the following implicit 
numerical scheme in these regions. Let introduce the grid Integration of equations (A3)-(A5) on ^ from ^, to ^,+1 results in the 
following expressions for density p„ velocity and the gas pressure Pg,, upstream of the forward shock > 1) at the grid knot 
with a number i in terms of p,+i and m,+i : 



,,/-'_ „old,,old 
Uifji - Pi Ui 



Rf 



uold 



Pi+lUi+i- 

-1 



(B.l) 



(B.2) 



rjold 



{ uold \^ 



R 



f ) 



-P 



- (jg - - Ui)T 



Pc,i+\ Pc,i 



1-3t 



(B.3) 



Here the superscript 'old' is referred to quantities at the previous time step and r is the step in the dimensionless time t. The 
quantity G; is given by the expression 



Gi=p1 



old 



f 



R 



f ) 



■ Pi+i - 



(B.4) 



In spite of impUcity Eqs (BI-B3) are readily solved recurrently from i = i^ax down to j = / corresponding to the position of the 
forward shock. The value of the forward shock radius Rf is related with the forward shock radius R°J^ at the previous time step as 
/?/ =/?^''^expTinEqs(BI-B4). 

Performing the integration of equations (A3)-(A5) on ^ from to ^, upstream of the reverse shock we obtain the expressions 
for for density, velocity p„ m, and and the gas pressure P^,, upstream of the reverse shock < I) in terms of p,_i and m,_i : 



^old 



UiH'. = uf + Ui- 



RbVfi^i-^i-i) 



RbVf(^i-^i-i) 



(B.5) 



Pi = Hi 



I + ■ ■ k 



1) 



8,1 



r)Old 



+ Po i-\ r 5 {y„ — I)(W,- - M;)t 

UJ ' RbVf{.ei-ei-,) ' " RbVf{^i-^i.,) 



(B.6) 



I + 3RfT 



0iui - Vb^i) + (jg - mfui - 0_,ui.i) 



The quantities Hi and H'^ are given by expressions 



(B.7) 



Hi=p1 



old 



+ P'-i TTTTTT^ — —s '^i = 1 + 



RbVfi^f-^U) 



RbVfi^i-^i-,) 



(B.8) 



Eqs (B5-B7) are solved recurrently from i = I to i = b corresponding to the position of the reverse shock. The value of the 
new reverse shock radius Rb in Eqs (B5-B8) was extrapolated using the reverse shock radius R'^'^ and velocity V^'^ as Rb = 
jfou _^_ v°'''Rj:T/Vf. 
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Appendix C. Solution of hydrodynamic equations in the downstream regions 

Hydrodynamical quantities are determined at the centers of the cells with the indexes / + 1 /2 in the downstream regions. We used 
an explicit numerical scheme of Trac & Pen ifisll for solution of Eqs. (A7)-(A9). In accordance with these equations we use the 
following variables 

1 2 2 1 3 2 

Ui+\I2 ^ Pi+l/2ri+i/2^f,b, f/i+1/2 = f/i+i/2Mi+l/2, f/i+i/2 = '";+i/2^/,/'ei+l/2, (C.l) 

and fluxes 

1 2 f b 2 2 f b 

^/+l/2 = ■^P-+l/2'";+i/2«,-;i/2' Fi+l/2 = ■^'"/+l/2(Pi+l/2«i+l/2M,-;i/2 + ^'g,/+l/2), 

3 1 fb 

= ■j7^''/+l/2(^'+I/2«;+i/2 + «/+l/2^j,i+I/2)- (C.2) 

The variables Uf^[^2 introduced as 

K[/2-'^^Ui,,,±Fl,^„ 1^1,2,3 (C.3) 



Here c, is the so called freezing speed 

1/2 



The values U"™i at the next instant of time are given by 



U^i2-Ul,^2-^(FL-Fj), (C.4) 



Km = - —{Fi, - Ff) + Tr,-,i/2-^ \2P,j,mHb - r,-.i/2^^^^^^j , (C.5) 

jTnew,3 tt3 _J'tr3 i?3x ^2 ,,, ^/ Pc,i+\ ~ Pc,i (n f,\ 

K-J,b yf K-f,b 

Here kf^b is the grid step downstream of the forward and reverse shock respectively. Note that a uniform grid in the downstream 
region was used. 

The fluxes F'. at the grid knots between discontinuities are given by 

f'i = 0.5{Fl'' + FT''), / = 1, 2, 3, (C.7) 
where the fluxes Ff'' at the grid knots are given by 

FP' = f/^;/2 + 0.5L(f/^;^, - Ut\/„ Ut^/, - f/^^/2)' (C.8) 

F?' - -u;Ai2 + ^-^m'm - v-Ai2^ v:t3,2 - v-Ai2)^ ^ = 1, 2, 3. (C.9) 

Here the function L(a, b) is the nonlinear flux limiter. 

The numerical scheme is reduced to the first order Godunov's one for L(a, b) = in Eqs. (C8), (C9). Such a scheme is very 
dissipative. The simplest second order scheme corresponds to L(a, b) - a. However it is unstable. One should use more complex 
flux limiters for stability. The different types of corresponding nonlinear limiters are available. Trac & Pen ifTsIl used a so-called 
Van-Leer limiter: 

L{a,b)^{ P^'f't^J' . (C.IO) 



0, ab 

We used a more dissipative "minmod" limiter: 

{a, ab > 0, \a\ < \b\, 
b, ab>0, \a\ >\b\, . (C.ll) 
0, ab < 0. 

For calculations of fluxes Fp' near discontinuities we also use L(a, b) - oi L(a, b) - a. 
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The fluxes fj. and fj, just downstream of the forward and reverse shock at the knots with numbers i - f and i = bare calculated 
using hydrodynamical quantities just upstream of the forward and reverse shocks: 

P),b = ■^Pf,bRlb(-'^f,b - Vf,b), Fji, = —Rj,,(pf^bUf,biuf,b - Vf,b) + Pg,f,b), 

^Ib = ^/Ibi^f-'^^f'" - ^f'"^ + '^f.i'PgXb)- (C. 12) 

The fluxes F' at the position of the contact discontinuity at the knot with number j = c are given by 
Fl = 0, = P'RlRf/Vf, Fl = V,Fl, 
where the pressure P' is found from the approximate solution of the Riemann problem for decay of an arbitrary discontinuity: 



, _ ^e+l/2 VPe-1/2 + ^c-1/2 s/p^+TJl + (Mc-1/2 - Me+1/2) yjjgPpc-lllPc+lll 

ypc-\ii + yjPc+i/i 

Here P = (Pc-i/i + Pc+i/2)/2. 

Hydrodynamical quantities at the new time step are found using Eq. (CI) after the calculation of new positions of the disconti- 
nuities (see Appendix D). 

For a stability of the numerical scheme considered the freezing speed Cj must be greater than the maximal sonic velocity. We use 
the foUowing expression for Cj 



Rf 

Cj = — max 



Q.5{Pc,i + Pc,M) + Pg,Mll , I _i 

\yg + l«'+l/2l ^fl 

I Pi+1/2 



,Rb<rM/2<Rf (C.14) 
The time step r was found from the relation r = 0.5 max(fc/, kb)/Cs. 



Appendix D. Calculation of the velocities of discontinuities 

The velocities of discontinuities are found from the approximate solutions of the Riemann problem describing a decay of an 
arbitrary discontinuity. For the forward and reverse shocks the solution for the shock velocities V/ and Vb are the following: 



Vf = Uf +p-/I^^L±^Pf_,i^ + ^Pf. (D.l) 



Vb = Ub-pl"^^ ^^Pb^i/2 + ^^Pb- (D.2) 

The shock velocities are determined by upstream values Pf,Pb,Uf,Ub,Pu,Pb and downstream values Pf-i/2,Pb+i/2 of the gas 
pressure, velocity and density. 

The expression for the velocity of the contact discontinuity Vc is given by 

Wc-l/2 VPc-1/2 + Mc+1/2 VPc+1/2 + '''^^'^J^^'^ - 

Vc = , ^ , — (D.3) 

VPc-l/2 + VPc+1/2 

Here P = (Pc-1/2 + Pc-,m)/2. 

The new values for the radii Rb, Rc,Rf = R°j^ exp(T) and real time t are now found: 



-old 



Rb = Rf + Rfr ^ 

* ^ Vf + Vf 



exp(0.5T), (D.4) 



y , yold 

Rc = Rf + Rfr— ^ exp(0.5T), (D.5) 



IRfr 

old . J_ 



t = f'^ + — - exp(0.5T). (D.6) 

y^ + yold 
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Appendix E. Solution of CR transport equation 

The finite difference scheme for Eqs (A6) and (A 10) can be written in the following form: 

Here are the values of n at the new time step calculated in the grid knot ^, or 77, and at (j. We use the following coefficients 
upstream of the shocks 



H = - - -J^^^lvi - ^ll/l) - I 4^ J. > 0. ' = ^'>l/2^^-l/2 - ^-l/2^i-U2, (E.4) 

_ _;^..3 _ .3 . /.ow _ ^ f -«r'' ^ 0' 5) 

The corresponding coefficients downstream of the shocks are given by 
2 

J _ r,j ""'-1/2 , f,h2 i 0, 77, > 0, 

/ _ ^- ^'+1/2 2 / 1, > 0, 

uJ j J r3 3 ^ di ( 4-h-\ di^O, ,22 ,T,o^ 

= - - ^^'■;+l/2 - ^i-l/l) - y I 4, J,. > 0. ' = 'T+l/2W!+l/2 - ^,-l/2Wi-l/2, (E.8) 

The forward and reverse shocks are situated at knots with numbers i = f and i = b respectively. The coefficients at these knots 
are given by 

^2 ^2 

4 = D{ , n-—^^ r, 4 = , ,,Rf ^'^^^'\ , (E.10) 



nRf.bJ^ - (w/ - Vf,i,mlb^f I i' f ! 1' , (E.3) 



Vf 4 — 

^/ = ~ 4 " 3^*^^/^/+l/2 - '■/-I/2) ^ ^/ = 4^/+l/2W/+l/2 - '•/-I/2W/-I/2, 



(E.ll) 



/ 



(R'ffuu2 - rlx,2>f - r/(Vf - ^f)Pf^6{j, jf) + -^nf , (E.12) 



and 

p2 2 



Vf 4 - 

K = - 4 - ^('■^+1/2 - ^^^i-1/2) ^^6. db = '-Ll/2VV'fc+l/2 - i^i^tl/2>»'i-l/2> (E.14) 



si = - RUlmK" - rfiwb - Vb)p/-^6U, k) + |nr^ . (E.15) 
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Here 6(j, is the Kronecker's symbol while jf and jh are numbers of the momentum grid knots corresponding to injection momenta 
Pf and ph respectively. 

The contact discontinuity is situated at the knot / - c. The coefficients at this knot are given by 

2 2 

J = ''''^'^ C-' = CE 16^ 



bi = -ai - ci - ^^(r^^j^2 - '"r-1/2) - J I 4^ >'o ' = f^+m^c^m - r,_^i2Wc-iii, (E.17) 



- ~^irT('"e+l/2 ~ '";-l/2)"c" j+\,old _^ J,old j ■ (E.18) 



37?/T^''^+i/2 '•,-i/2^«c 3/4 «f -wf''', ^/c.>0 

We use a standard method to solve the tridiagonal set of equations (El) (e.g. Godunov Jl^). At first the following coefficients 
L^^j^2 ^nd ^/|_i^2 ^'^^ recurrently calculated from / = 1 up to / = i,„ax — 1: 



g'. - a'.K'. , 

n = = '. . ' (E.19) 



/2 



i ;-l/2 / t t-i/2 I 

The values Lj^^ = 1 and ^Tj^^ = correspond to the boundary condition - at the center of the remnant. The functions are 
calculated recurrently from / - i^ax - 1 down to ; = as 

«/ = ^L/2<i+<i/2- (E.20) 

We put «^ = at / = imax- This describes an absorbtion of particles at the boundary of the simulation domain. The equations (El 9), 
(E20) are recurrently used for all values of j starting with the lowest energies at j - jmin + 1 ■ 



Appendix F. Initial conditions and the code performance 

We use a uniform spacial grid downstream of the shocks and the following non-uniform grids upstream of the forward shock 

= 1 - ^min + ^min exp[(; - f)ki], ki = \n[(^„,ax + ^min ' ^) I ^min] I (imax ' f) (F-l) 

and upstream of the reverse shock 

= 2/^4 + [in (1 + (2^„„„)('-2*4)) _ + 2^„„„)] / ln(2^„„„), = 1 jb. (F.2) 

Here the parameter ^„,ax describes the ratio of radii of the simulation domain and forward shock. This parameter is constant during 
the simulation. We use the value ^„,ax - 2. The value of the parameter ^„„>, ~ 10"" should be small enough to resolve the spacial 
variation scale of the lowest energy cosmic rays upstream of the shocks. 

We use the following profiles of the gas density, velocity and pressure at the initial moment of time f = 0: 

p — Pq, u - 0, Pg - Pq, r > Rf, 

p = 4po, u = 0J5Vf, Pg - OJSpoVj, Rb<r< Rf, 

p = po{rlRb)-\ u = \.5VfrlRb, Pg = 0, r,j < r < Rb, *■ ' 

p = po{rejlRb)-\ u = l.5Vfr/Rb,Pg = 0, r < r,j. 

Here po ™d f are the gas density and pressure in the circumstellar medium while r^y - IVejRb/^Vf is the initial radius of the 
flat part of the ejecta density distribution. The initial velocities of the reverse shock and contact discontinuity are Vb = O.SV/ and 
Vc = Q.lSVf respectively. Their initial radii can be chosen rather arbitrary. We use values Rb = 0.9/?/ and R^ - 0.95Rf. 

The initial cosmic ray number density is zero. At every time step the system of equations (l)-(4) is solved in the following order: 

1) The hydrodynamical equations are solved in the upstream regions (see Appendix B). 

2) These equations are solved in the downstream regions and the speeds of discontinuities are determined (see Appendix C and 
D). 

3) The parameters of CR transport depending on the hydrodynamical quantities are calculated. 

4) CR transport equation is solved (see appendix E). Cosmic ray pressure is calculated. 

It is important to note that we use the finite difference method for solution of CR transport equations even at the shock positions 
at grid knots / = / and ; = b (see Eqs. (E10-E15)). These equations contain terms originating from the first derivative on time (see 
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terms containing r in Eqs.(E10-E15)). These terms will not appear if we use exact boundary conditions at the shock fronts which 
can be derived by integration of Eq. (4) in the vicinity of the shocks. Our method permits to use relatively large time steps r while 
the use of the exact boundary conditions results in a numerical instability for the large time steps. Physically this is because the low 
energy particles have very small acceleration times in the case of Bohm-like diffusion. This results in the numerical instabilities. It 
is possible to use many iterations at every time step 101 or to use very small time steps near the shock front ^ to avoid this problem. 
This results in significantly longer times of simulations. 

The numerical scheme with Eqs. (E10-E15) artificially increases acceleration time of low-energy particles and makes the calcu- 
lations to be stable. This means that our method can not be used for simulations of fast varying processes. However if the shocks 
propagate in the smooth environment then the low-energy particles are accelerated in a quasi-steady regime and our method is well 
justified. 

The numerical results shown in Figs (I)-(IO) were obtained using the grid with 200x200 radius-momentum cells in the every 
of four regions (two upstream and two downstream regions of the shocks). The initial radius of the flat part of the ejecta density 
distribution r^j is r^j = Rh/(!>- The calculation of the remnant evolution takes two hours at PC. The total energy is conserved with 
5% accuracy. Preliminary results can be obtained using ten times faster crude calculations with 100x100 grid. 
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